#=============================================================================#
#
# PROJECT:        Does Islamist Terrorism Still Affect Political Attitudes?
# AUTHORS:        ** anonymized for review **
# CONTACT:        ** anonymized for review **
# LAST MODIFIED:  November 28, 2025
# 
#=============================================================================#
#
# This R file contains the code used to replicate the robustness analyses
# 
#=============================================================================#

rm(list=ls())
getwd()

# Install and load all necessary packages with the ipak function:
# i.e., check to see if packages are installed. 
# Install them if they are not, 
# then load them into the R session.

ipak <- function(pkg){  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
sapply(pkg, require, character.only = TRUE)
}

packages <- c("readxl", "metafor", "metaSEM", "ggplot2",
              "tidyverse", "dplyr", "splines",
              "modelsummary", "stargazer", "texreg") 
ipak(packages)


# Data ------------------------------------------------------------------------

## Data: Full dataset, now including attacks outside the west
df <- readRDS("data/df.rds") 
df <- df %>%
  mutate(StudyYear = StudyYear - 2001)


# Robustness checks -----------------------------------------------------------

## Robustness checks of Model 2 of over-time variation

# Check 1: Including Studies on Non-Western Attacks
m3.r1 <- rma.mv(d, var_d,
                mods = ~ StudyYear + # StudyYear as moderator
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample +
                  Outgroup + Rally +
                  US,
             random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
             data = df)
m3.r1

# now re-delete 3 studies on effects of terrorism conducted outside Western democracies 
## for other robustness checks (to be consistent with sample used in main models)
df <- df %>% 
  filter((ExactAttack_2 != "2002 Bali bombings" &  # attacks on western countries 
            ExactAttack_2 != "2008 Mumbai attacks" &  # attacks on western countries 
            ExactAttack_2 != "Gaza - Israel conflict") %>%
           replace_na(TRUE))


# Check 2: Including 9/11 Dummy
m3.r2 <- rma.mv(d, var_d,
                mods = ~ StudyYear + # StudyYear as moderator
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample +
                  Outgroup + Rally +
                  US + NineEleven,
                random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                data = df)
m3.r2


# Check 3: Excluding Trait-Like Predispositions
m3.r3 <- rma.mv(d, var_d,
                mods = ~ StudyYear + # StudyYear as moderator
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample +
                  Outgroup + Rally +
                  US,
                random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                data = df[df$PA_Category != 5 & # RWA
                            df$PA_Category != 6 & # SDO
                            df$PA_Category != 7 & # General ideology
                            df$PA_Category != 11,]) # Mixed or unclear
m3.r3


# Check 4: Restricting to Natural Experiments
m3.r4 <- rma.mv(d, var_d,
                mods = ~ StudyYear + # StudyYear as moderator
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample +
                  Outgroup + Rally +
                  US,
                random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                data = df[df$TypeStudy == 2,]) # Natural experiments only
m3.r4


# Check 5: Excluding Multi-Variate Regression Estimates
m3.r5 <- rma.mv(d, var_d,
                mods = ~ StudyYear + # StudyYear as moderator
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample +
                  Outgroup + Rally +
                  US,
                random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                data = df[df$Regression == 0,]) # No regression-based estimates
m3.r5


# Check 6: Excluding Outliers
# define outliers
df <- df %>% 
  mutate(Outlier = ifelse(abs(d - mean(d)) > 3*sd(d), 1, 0))
# exclude outliers
m3.r6 <- rma.mv(d, var_d,
                mods = ~ StudyYear + # StudyYear as moderator
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample +
                  Outgroup + Rally +
                  US,
                random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                data = df[df$Outlier == 0,]) # No outliers
m3.r6


# Table S1.3

# Step 1: Create a named list of models
models <- list(
  "(1)" = m3.r1,
  "(2)" = m3.r2,
  "(3)" = m3.r3,
  "(4)" = m3.r4,
  "(5)" = m3.r5,
  "(6)" = m3.r6
)

# Create the table
cm <- c("StudyYear"     = "Year of study",
        "CasualtiesLess101"  = "Fatalities: $<$10",
        "Casualties10to1001" = "Fatalities: 10–100",
        "Casualties100plus1" = "Fatalities: $>$100",
        "Exp1"          = "Randomized experiment",
        "NatExp1"       = "Natural experiment",
        "ProbSample1"   = "Probability sample",
        "Outgroup"      = "Outgroup hostility",
        "Rally"         = "Rally tendencies",
        "US1"           = "US dummy",
        "NineEleven1"   = "9/11 dummy",
        "intercept"     = "Constant")
modelsummary(models, 
             coef_map = cm,
             stars = c("+" = .1, "*" = .05, "**" = .01, "***" = .001),
             fmt = 3,
             gof_omit = 'IC|Log|Adj|AIC|BIC|RMSE',
             escape = FALSE,
             output = "latex_tabular")



## Robustness check of past attacks -------------------------------------------

## Get estimates for...
## Get estimates for...
# ... Model 2: Attacks in the same country during the 5 years prior to the study.
m2_past_1 <- rma.mv(d, var_d, 
                    mods = ~ log1p(PastAttacks_Country_5yr) + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_1

# ... Model 2: Attacks in the West during the 5 years prior to the study.
m2_past_2 <- rma.mv(d, var_d, 
                    mods = ~ log1p(PastAttacks_West_5yr) + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_2

# ... Model 2: Attacks in the same country, no time restriction.
m2_past_3 <- rma.mv(d, var_d, 
                    mods = ~ log1p(PastAttacks_Country_All) + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_3

# ... Model 2: Attacks in the West, no time restriction.
m2_past_4 <- rma.mv(d, var_d, 
                    mods = ~ log1p(PastAttacks_West_All) + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      Outgroup + Rally +
                      US,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_4


## Table S1.4

# Step 1: Create a named list of models
models <- list(
  "Model 1" = m2_past_1,
  "Model 2" = m2_past_2,
  "Model 2" = m2_past_3,
  "Model 4" = m2_past_4
)

# Create the table
cm <- c("log(PastAttacks_Country_5yr + 1)" = "Log(Number of previous Islamist terrorist attacks (same country, 5 years))+1",
        "log(PastAttacks_West_5yr + 1)" = "Log(Number of previous Islamist terrorist attacks (all Western countries, 5 years))+1",
        "log(PastAttacks_Country_All + 1)" = "Log(Number of previous Islamist terrorist attacks (same country, all years))+1",
        "log(PastAttacks_West_All + 1)" = "Log(Number of previous Islamist terrorist attacks (all Western countries, all years))+1",
        "intercept"     = "Constant")
modelsummary(models,
             coef_map = cm,
             stars = c("+" = .1, "*" = .05, "**" = .01, "***" = .001),
             fmt = 3,
             gof_omit = 'IC|Log|Adj|AIC|BIC|RMSE',
             escape = FALSE,
             output = "latex_tabular")



# Complementary Analyses ------------------------------------------------------

## Outcome-specific models ----------------------------------------------------

## Get estimates for...
# ... Sub-sample: Conservative Shift
# ... Model 2
m3.con <- rma.mv(d, var_d, 
                 mods = ~ StudyYear + # StudyYear as moderator
                   Exp + NatExp +
                   ProbSample + 
                   CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                   US,
                 random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
             data = df[df$Conservative==1,])
m3.con


## Get estimates for...
# ... Sub-sample: Outgroup Hostility
# ... Model 2
m3.out <- rma.mv(d, var_d, 
                 mods = ~ StudyYear + # StudyYear as moderator
                   Exp + NatExp +
                   ProbSample + 
                   CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                   US,
                 random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                 data = df[df$Outgroup==1,])
m3.out


## Get estimates for...
# ... Sub-sample: Rally/Activation
# ... Model 2
m3.rally <- rma.mv(d, var_d, 
                 mods = ~ StudyYear + # StudyYear as moderator
                   Exp + NatExp +
                   ProbSample + 
                   CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                   US,
                 random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                 data = df[df$Rally==1,])
m3.rally


## Table S1.4

# Step 1: Create a named list of models
models <- list(
  "Conservatism" = m3.con,
  "Outgroup" = m3.out,
  "Rally" = m3.rally
)

# Create the table
cm <- c("StudyYear"     = "Year of study",
        "CasualtiesLess101"  = "Fatalities: $<$10",
        "Casualties10to1001" = "Fatalities: 10–100",
        "Casualties100plus1" = "Fatalities: $>$100",
        "Exp1"          = "Randomized experiment",
        "NatExp1"       = "Natural experiment",
        "Rally"         = "Rally tendencies",
        "US1"           = "US dummy",
        "intercept"     = "Constant")
modelsummary(models, 
             coef_map = cm,
             stars = c("+" = .1, "*" = .05, "**" = .01, "***" = .001),
             fmt = 3,
             gof_omit = 'IC|Log|Adj|AIC|BIC|RMSE',
             escape = FALSE,
             output = "latex_tabular")



## Non-Linear Dynamics --------------------------------------------------------

#### --- Non-linear over-time models ----------------------

# 1) Settings ---
ctrls <- c("CasualtiesLess10","Casualties10to100","Casualties100plus",
           "Exp","NatExp","ProbSample","Outgroup","Rally","US")

# recode all binary 0/1 factors to numeric
df[ctrls] <- lapply(df[ctrls], function(x) as.numeric(as.character(x)))

# compute control means (assumes numeric 0/1; adjust here if you have factors)
ctrl_means <- sapply(df[ctrls], function(x) mean(x, na.rm = TRUE))

# 2) Fit four over-time models ---
m_year_lin <- rma.mv(d, var_d,
                     mods = ~ StudyYear + 
                       CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                       Exp + NatExp +
                       ProbSample +
                       Outgroup + Rally +
                       US,
                     random = ~ 1 | ID_R/ID_ES_Unique, data = df)

m_year_quad <- rma.mv(d, var_d,
                      mods = ~ StudyYear + I(StudyYear^2) +
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp +
                        ProbSample +
                        Outgroup + Rally +
                        US,
                      random = ~ 1 | ID_R/ID_ES_Unique, data = df)

m_year_cubic <- rma.mv(d, var_d,
                       mods = ~ StudyYear + I(StudyYear^2) + I(StudyYear^3) +
                         CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                         Exp + NatExp +
                         ProbSample +
                         Outgroup + Rally +
                         US,
                       random = ~ 1 | ID_R/ID_ES_Unique, data = df)

m_year_spline <- rma.mv(d, var_d,
                        mods = ~ ns(StudyYear, df = 3) + 
                          CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                          Exp + NatExp +
                          ProbSample +
                          Outgroup + Rally +
                          US,
                        random = ~ 1 | ID_R/ID_ES_Unique, data = df)


# coefficient table for Linear, Quadratic, Cubic, Spline
texreg::texreg(
  list(m_year_lin, m_year_quad, m_year_cubic, m_year_spline),
  custom.model.names = c("Linear","Quadratic","Cubic","Spline (df=3)"),
  caption  = "Over-time meta-regressions with non-linear specifications",
  label    = "tab:overtime_nonlin_models",
  dcolumn  = TRUE,
  use.packages = FALSE,    # keeps table self-contained (no extra \usepackage)
  float.pos = "H",
  file = "output/tables/overtime_nonlin_models.tex"
)

# 3) Prediction grid (controls at means ---
year_seq <- seq(min(df$StudyYear, na.rm = TRUE),
                max(df$StudyYear, na.rm = TRUE), by = 1)

new_year <- data.frame(StudyYear = year_seq)
for (v in ctrls) new_year[[v]] <- ctrl_means[[v]]
new_year

# model matrices must align with each fitted model (drop intercept column)
X_lin <- model.matrix(~ StudyYear + 
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp + ProbSample + Outgroup + Rally + US,
                      data = new_year)[, -1]

X_quad <- model.matrix(~ StudyYear + I(StudyYear^2) + 
                         CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                         Exp + NatExp + ProbSample + Outgroup + Rally + US,
                       data = new_year)[, -1]

X_cubic <- model.matrix(~ StudyYear + I(StudyYear^2) + I(StudyYear^3) +
                          CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                          Exp + NatExp + ProbSample + Outgroup + Rally + US,
                        data = new_year)[, -1]

X_spline <- model.matrix(~ ns(StudyYear, df = 3) +
                           CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                           Exp + NatExp + ProbSample + Outgroup + Rally + US,
                         data = new_year)[, -1]

# 4) Predict and bind into one DF ---
p_lin    <- predict(m_year_lin,    newmods = X_lin)
p_quad   <- predict(m_year_quad,   newmods = X_quad)
p_cubic  <- predict(m_year_cubic,  newmods = X_cubic)
p_spline <- predict(m_year_spline, newmods = X_spline)

pred_year <- dplyr::bind_rows(
  tibble::tibble(Model = "Linear",   StudyYear = new_year$StudyYear,
                 pred = p_lin$pred,  ci.lb = p_lin$ci.lb,  ci.ub = p_lin$ci.ub),
  tibble::tibble(Model = "Quadratic",StudyYear = new_year$StudyYear,
                 pred = p_quad$pred, ci.lb = p_quad$ci.lb, ci.ub = p_quad$ci.ub),
  tibble::tibble(Model = "Cubic",    StudyYear = new_year$StudyYear,
                 pred = p_cubic$pred,ci.lb = p_cubic$ci.lb,ci.ub = p_cubic$ci.ub),
  tibble::tibble(Model = "Spline (df=3)", StudyYear = new_year$StudyYear,
                 pred = p_spline$pred,ci.lb = p_spline$ci.lb,ci.ub = p_spline$ci.ub)
)

# 5) Plot ---
p <- ggplot(pred_year, aes(StudyYear, pred)) +
  geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = 0.15) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ factor(Model, c("Linear", "Quadratic", "Cubic", "Spline (df=3)")), 
             ncol = 2) +
  labs(x = "Study Year (centered; 0 = 2001)",
       y = "Predicted effect size (Cohen's d)") +
  theme_bw(base_size = 12) +
  coord_cartesian(ylim = c(-0.1, 0.4)) 
print(p)
ggplot2::ggsave("output/figures/meta_nonlinear.pdf", p, width = 7.0, height = 5.0)


# 6) Quick fit comparison printed to console ---
## refit ML versions for valid LRTs (same random structure)
m_year_lin_ml   <- update(m_year_lin,   method = "ML")
m_year_quad_ml  <- update(m_year_quad,  method = "ML")
m_year_cubic_ml <- update(m_year_cubic, method = "ML")
m_year_spl_ml   <- update(m_year_spline,method = "ML")  # for AIC/BIC only (not nested)


#### --- Non-linear attack models ----------------------

# 1) Transformations ---
df <- df %>%
  mutate(
    logPast_W  = log1p(PastAttacks_Country_All),
    sqrtPast_W = sqrt(PastAttacks_Country_All)
  )

# 2) Fit models (REML for coefficients/plots) ---
m_att_lin   <- rma.mv(d, var_d,
                      mods = ~ PastAttacks_Country_All + 
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp + ProbSample + Outgroup + Rally + US,
                      random = ~ 1 | ID_R/ID_ES_Unique, data = df, method = "REML")

m_att_quad  <- rma.mv(d, var_d,
                      mods = ~ PastAttacks_Country_All + I(PastAttacks_Country_All^2) + 
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp + ProbSample + Outgroup + Rally + US,
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus + US,
                      random = ~ 1 | ID_R/ID_ES_Unique, data = df, method = "REML")

m_att_log   <- rma.mv(d, var_d,
                      mods = ~ logPast_W + 
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp + ProbSample + Outgroup + Rally + US,
                      random = ~ 1 | ID_R/ID_ES_Unique, data = df, method = "REML")

m_att_sqrt  <- rma.mv(d, var_d,
                      mods = ~ sqrtPast_W + 
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp + ProbSample + Outgroup + Rally + US,
                      random = ~ 1 | ID_R/ID_ES_Unique, data = df, method = "REML")

m_att_spl   <- rma.mv(d, var_d,
                      mods = ~ ns(PastAttacks_Country_All, df = 3) + 
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                        Exp + NatExp + ProbSample + Outgroup + Rally + US,
                      random = ~ 1 | ID_R/ID_ES_Unique, data = df, method = "REML")

# 3) Coefficient table ---
texreg::texreg(list(m_att_lin, m_att_quad, m_att_log, m_att_sqrt, m_att_spl),
               custom.model.names = c("Linear","Quadratic","Log","Sqrt","Spline (df=3)"),
               caption = "Exposure (PastAttacks\\_West\\_All) meta-regressions with non-linear specifications.",
               label = "tab:attacks_nonlinear",
               dcolumn = TRUE, use.packages = FALSE, float.pos = "H",
               file = "output/tables/exposure_nonlin_models.tex")

# 4) Prediction grid ---
att_seq <- seq(min(df$PastAttacks_Country_All, na.rm = TRUE),
               max(df$PastAttacks_Country_All, na.rm = TRUE), length.out = 150)

new_att <- data.frame(
  PastAttacks_Country_All = att_seq,
  logPast_W  = log1p(att_seq),
  sqrtPast_W = sqrt(att_seq)
)

# hold controls at their means
for (v in ctrls) new_att[[v]] <- ctrl_means[[v]]

# Build design matrices and predict for each model (drop intercept)
X_lin <- model.matrix(~ PastAttacks_Country_All + Exp + NatExp + ProbSample + Outgroup + Rally +
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus + US,
                      data = new_att)[, -1]
X_quad <- model.matrix(~ PastAttacks_Country_All + I(PastAttacks_Country_All^2) + Exp + NatExp + ProbSample +
                         Outgroup + Rally + CasualtiesLess10 + Casualties10to100 +
                         Casualties100plus + US,
                       data = new_att)[, -1]
X_log <- model.matrix(~ logPast_W + Exp + NatExp + ProbSample + Outgroup + Rally +
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus + US,
                      data = new_att)[, -1]
X_sqrt <- model.matrix(~ sqrtPast_W + Exp + NatExp + ProbSample + Outgroup + Rally +
                         CasualtiesLess10 + Casualties10to100 + Casualties100plus + US,
                       data = new_att)[, -1]
X_spl <- model.matrix(~ ns(PastAttacks_Country_All, df = 3) + Exp + NatExp + ProbSample + Outgroup + Rally +
                        CasualtiesLess10 + Casualties10to100 + Casualties100plus + US,
                      data = new_att)[, -1]

p_lin   <- predict(m_att_lin,   newmods = X_lin)
p_quad  <- predict(m_att_quad,  newmods = X_quad)
p_log   <- predict(m_att_log,   newmods = X_log)
p_sqrt  <- predict(m_att_sqrt,  newmods = X_sqrt)
p_spl   <- predict(m_att_spl,   newmods = X_spl)

pred_att <- bind_rows(
  tibble(Model = "Linear",  PastAttacks = new_att$PastAttacks_Country_All,
         pred = p_lin$pred, ci.lb = p_lin$ci.lb, ci.ub = p_lin$ci.ub),
  tibble(Model = "Quadratic", PastAttacks = new_att$PastAttacks_Country_All,
         pred = p_quad$pred, ci.lb = p_quad$ci.lb, ci.ub = p_quad$ci.ub),
  tibble(Model = "Log", PastAttacks = new_att$PastAttacks_Country_All,
         pred = p_log$pred, ci.lb = p_log$ci.lb, ci.ub = p_log$ci.ub),
  tibble(Model = "Sqrt", PastAttacks = new_att$PastAttacks_Country_All,
         pred = p_sqrt$pred, ci.lb = p_sqrt$ci.lb, ci.ub = p_sqrt$ci.ub),
  tibble(Model = "Spline (df=3)", PastAttacks = new_att$PastAttacks_Country_All,
         pred = p_spl$pred, ci.lb = p_spl$ci.lb, ci.ub = p_spl$ci.ub)
)

# 5) Plot ---
p_exp <- ggplot(pred_att, aes(PastAttacks, pred)) +
  geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), alpha = 0.15) +
  geom_line(linewidth = 0.9) +
  facet_wrap(~ factor(Model, c("Linear", "Quadratic",  "Log", "Sqrt", "Spline (df=3)")), 
                      ncol = 2) +
  labs(x = "Number of past attacks in country of study",
       y = "Predicted effect size (Cohen's d)") +
  theme_bw(base_size = 12) +
  coord_cartesian(ylim = c(-0.1, 0.4)) 
print(p_exp)

ggplot2::ggsave("output/figures/attacks_nonlinear.pdf", p_exp, width = 7, height = 6)


## Proportion of null results over time ------------------------------------

df$NullResult <- as.factor(ifelse(df$pval_d > 0.05, "Insignificant result", "Significant result"))
summary(df$NullResult)

# adjust the margin settings
par(mar = c(2, 2, 0, 0))

# create the mosaic plot
mosaicplot(StudyYear ~ NullResult, data = df, 
           col = c('gray80', 'gray20'),
           ylab = "Proportions",
           xlab = "Year of data collection",
           cex.axis = 0.9, las = 3,
           main = NULL)

# reset the margin settings to default
par(mar = c(5, 4, 4, 2) + 0.1)

# save the plot as a PDF file
meta.mosaic <- recordPlot()
pdf(file="output/figures/meta_null_results.pdf",
    width = 7.5, height = 5.5)
meta.mosaic
dev.off()



## Heterogeneity ------------------------------------
m2_mods <- rma.mv(d, var_d, 
                mods = ~ StudyYear*FemaleP + 
                  StudyYear*AgeMean + 
                  StudyYear*StudentPop +
                  StudyYear*US + 
                  CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                  Exp + NatExp +
                  ProbSample + StudentPop +
                  Outgroup + Rally +
                  US +
                  FemaleP +
                  AgeMean,
                random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                data = df)
m2_mods


m2_past_1_mods <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_Country_5yr*FemaleP + 
                      PastAttacks_Country_5yr*AgeMean + 
                      PastAttacks_Country_5yr*StudentPop +
                      PastAttacks_Country_5yr*US + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      StudentPop +
                      Outgroup + Rally +
                      US +
                      FemaleP +
                      AgeMean,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_1_mods

# ... Model 2: Attacks in the West during the 5 years prior to the study.
m2_past_2_mods <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_West_5yr*FemaleP + 
                      PastAttacks_West_5yr*AgeMean + 
                      PastAttacks_West_5yr*StudentPop +
                      PastAttacks_West_5yr*US + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      StudentPop +
                      Outgroup + Rally +
                      US +
                      FemaleP +
                      AgeMean,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_2_mods

# ... Model 2: Attacks in the same country, no time restriction.
m2_past_3_mods <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_Country_All*FemaleP + 
                      PastAttacks_Country_All*AgeMean + 
                      PastAttacks_Country_All*StudentPop +
                      PastAttacks_Country_All*US + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      StudentPop +
                      Outgroup + Rally +
                      US +
                      FemaleP +
                      AgeMean,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_3_mods

# ... Model 2: Attacks in the West, no time restriction.
m2_past_4_mods <- rma.mv(d, var_d, 
                    mods = ~ PastAttacks_West_All*FemaleP + 
                      PastAttacks_West_All*AgeMean + 
                      PastAttacks_West_All*StudentPop +
                      PastAttacks_West_All*US + 
                      CasualtiesLess10 + Casualties10to100 + Casualties100plus +
                      Exp + NatExp +
                      ProbSample +
                      StudentPop +
                      Outgroup + Rally +
                      US +
                      FemaleP +
                      AgeMean,
                    random = ~ 1 | ID_R/ID_ES_Unique, # three-level model
                    data = df)
m2_past_4_mods


# Table S1.X

# Step 1: Create a named list of models
models <- list(
  "(1)" = m2_mods,
  "(2)" = m2_past_1_mods,
  "(3)" = m2_past_2_mods,
  "(4)" = m2_past_3_mods,
  "(5)" = m2_past_4_mods
)

# Create the table
cm <- c("intrcpt"       = "Constant",
        "StudyYear"     = "Year of study",
        "PastAttacks_Country_5yr" = "Same country in previous 5 years",
        "PastAttacks_West_5yr" = "All Western countries in previous 5 years",
        "PastAttacks_Country_All" = "Same country in all previous years",
        "PastAttacks_West_All" = "All Western countries in all previous years",
        "FemaleP"       = "Percentage of Female Respondents",
        "AgeMean"       = "Mean Age of Respondents",
        "CasualtiesLess101"  = "Fatalities: $<$10",
        "Casualties10to1001" = "Fatalities: 10–100",
        "Casualties100plus1" = "Fatalities: $>$100",
        "Exp1"          = "Randomized experiment",
        "NatExp1"       = "Natural experiment",
        "ProbSample1"   = "Probability sample",
        "StudentPop"    = "Student sample",
        "Outgroup"      = "Outgroup hostility",
        "Rally"         = "Rally tendencies",
        "US"            = "US study",
        
        "StudyYear:FemaleP" = "Year of Study * Female share",
        "StudyYear:AgeMean" = "Year of Study  * Mean age",
        "StudyYear:StudentPop" = "Year of Study  * Student population",
        "StudyYear:US" = "Year of Study  * US study",
        
        "PastAttacks_Country_5yr:FemaleP" = "Past attacks * Female share",
        "PastAttacks_Country_5yr:AgeMean" = "Past attacks * Mean age",
        "PastAttacks_Country_5yr:StudentPop" = "Past attacks * Student population",
        "PastAttacks_Country_5yr:US" = "Past attacks *  US study",
        
        "PastAttacks_West_5yr:FemaleP" = "Past attacks * Female share",
        "PastAttacks_West_5yr:AgeMean" = "Past attacks * Mean age",
        "PastAttacks_West_5yr:StudentPop" = "Past attacks * Student population",
        "PastAttacks_West_5yr:US" = "Past attacks *  US study",
        
        "PastAttacks_Country_All:FemaleP" = "Past attacks * Female share",
        "PastAttacks_Country_All:AgeMean" = "Past attacks * Mean age",
        "PastAttacks_Country_All:StudentPop" = "Past attacks * Student population",
        "PastAttacks_Country_All:US" = "Past attacks *  US study",
        
        "PastAttacks_West_All:FemaleP" = "Past attacks * Female share",
        "PastAttacks_West_All:AgeMean" = "Past attacks * Mean age",
        "PastAttacks_West_All:StudentPop" = "Past attacks * Student population",
        "PastAttacks_West_All:US" = "Past attacks *  US study")
modelsummary(models, 
             coef_map = cm,
             stars = c("*" = .05, "**" = .01, "***" = .001),
             fmt = 3,
             gof_omit = 'IC|Log|Adj|AIC|BIC|RMSE',
             escape = FALSE,
             output = "latex_tabular")



